Non-Markovian dynamics of a qubit coupled to an Ising spin bath 
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We study the analytically solvable Ising model of a single qubit system coupled to a spin bath. The 
purpose of this study is to analyze and elucidate the performance of Markovian and non-Markovian 
master equations describing the dynamics of the system qubit, in comparison to the exact solution. 
We find that the time-convolutionless master equation performs particularly well up to fourth order 
in the system-bath coupling constant, in comparison to the Nakajima-Zwanzig master equation. 
Markovian approaches fare poorly due to the infinite bath correlation time in this model. A recently 
proposed post-Markovian master equation performs comparably to the time-convolutionless master 
equation for a properly chosen memory kernel, and outperforms all the approximation methods 
considered here at long times. Our findings shed light on the applicability of master equations to 
the description of reduced system dynamics in the presence of spin-baths. 
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I. INTRODUCTION 

A major conceptual as well as technical difficulty in the 
practical implementation of quantum information pro- 
cessing and quantum control schemes is the unavoidable 
interaction of quantum systems with their environment. 
This interaction can destroy quantum superpositions and 
lead to an irreversible loss of information, a process gener- 
ally known as decoherence. Understanding the dynamics 
of open quantum systems is therefore of considerable im- 
portance. The Schrodinger equation, which describes the 
evolution of closed systems, is generally inapplicable to 
open systems, unless one includes the environment in the 
description. This is, however, generally difficult, due to 
the large number of environment degrees of freedom. An 
alternative is to develop a description for the evolution 
of only the subsystem of interest. A multitude of dif- 
ferent approaches have been developed in this direction, 
exact as well as approximate [l|, [2] • Typically the exact 
approaches are of limited practical usefulness as they are 
either phenomenological or involve complicated integro- 
differential equations. The various approximations lead 
to regions of validity that have some overlap. Such tech- 
niques have been studied for many different models, but 
their performance in general, is not fully understood. 

In this work we consider an exactly solvable model of 
a single qubit (spin 1/2 particle) coupled to an environ- 
ment of qubits. We are motivated by the physical im- 
portance of such spin bath models @ in the description 
of decoherence in solid state quantum information pro- 
cessors, such as systems based on the nuclear spin of 
donors in semiconductors 0, Q , or on the electron spin 
in quantum dots Q. Rather than trying to accurately 
model decoherence due to the spin bath in such systems 
(as in, e.g., Refs. @, Q)) our g° a l m this work is to com- 
pare the performance of different master equations which 
have been proposed in the literature. Because the model 
we consider is exactly solvable, we are able to accurately 



assess the performance of the approximation techniques 
that we study. In particular, we study the Born-Markov 
and Born master equations, and the perturbation expan- 
sions of the Nakajima-Zwanzig (NZ) and the time- 
convolutionless (TCL) master equations [ll|, Q3 up to 
fourth order in the coupling constant. We also study the 
post-Markovian (PM) master equation proposed in [l3[ . 

The dynamics of the system qubit in the model we 
study is highly non-Markovian and hence we do not ex- 
pect the traditional Markovian master equations com- 
monly used, e.g., in quantum optics [14j and nuclear 
magnetic resonance [15|, to be accurate. This is typ- 
ical of spin-baths, and was noted, e.g., by Breuer et 
al. [l6j]. The work by Breuer et al. (as well as by 
other authors in a number of subsequent publications 
0, [H, EI HI El, HI) is conceptually close to ours 
in that in both cases an analytically solvable spin-bath 
model is considered and the analytical solution for the 
open system dynamics is compared to approximations. 
However, there are also important differences, namely, in 
Ref. [HI a so-called spin-star system was studied, where 
the system spin has equal couplings to all the bath spins, 
and these are of the XY exchange-type. In contrast, in 
our model the system spin interacts via Ising couplings 
with the bath spins, and we allow for arbitrary coupling 
constants. As a result there are also important differ- 
ences in the dynamics. For example, unlike the model in 
Ref. [lj| , for our model we find that the odd order terms 
in the perturbation expansions of Nakajima-Zwanzig and 
time-convolutionless master equations are non- vanishing. 
This reflects the fact that there is a coupling between the 
x and y components of the Bloch vector which is absent 
in (l6| . In view of the non-Markovian behavior of our 
model, we also discuss the relation between a represen- 
tation of the analytical solution of our model in terms of 
completely positive maps, and the Markovian limit ob- 
tained via a coarse-graining method introduced in [23j |. 
and the performance of the post-Markovian master equa- 
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tion [l3|. 

This paper is organized as follows. In Sec II, we present 
the model, derive the exact solution and discuss its be- 
havior in the limit of small times and large number of 
bath spins, and in the cases of discontinuous spectral den- 
sity co-domain and alternating sign of the system-bath 
coupling constants. In Sec. Ill, we consider second order 
approximation methods such as the Born-Markov and 
Born master equations, and a coarse-graining approach 
to the Markovian semigroup master equation. Then we 
derive solutions to higher order corrections obtained from 
the Nakajima-Zwanzig and time-convolutionless projec- 
tion techniques as well as derive the optimal approxima- 
tion achievable through the post-Markovian master equa- 
tion. In Sec. IV, we compare these solutions for various 
parameter values in the model and plot the results. Fi- 
nally in Sec. V, we present our conclusions. 



II. EXACT DYNAMICS 

A. The model 

We consider a single spin-i system (i.e., a qubit with 
a two-dimensional Hilbert space Hs) interacting with a 
bath of N spin-i particles (described by an ./V-fold tensor 
product of two-dimensional Hilbert spaces denoted Hb)- 
We model the interaction between the system qubit and 
the bath by the Ising Hamiltonian 

N 

H'j = aa z ®Y,9n<, (1) 
n=i 

where g n are dimensionless real-valued coupling con- 
stants in the interval [—1,1], and a > is a parame- 
ter having the dimension of frequency (we work in units 
in which h = 1), which describes the coupling strength 
and will be used below in conjunction with time (at) for 
perturbation expansions. The system and bath Hamilto- 
nians are 

H S = ^w a z (2) 

and 

N 1 

Hb = J2 o n "<- ( 3 ) 

n=l 

For definiteness, we restrict the frequencies u>o and fl n to 
the interval [—1,1], in inverse time units. Even though 
the units of time can be arbitrary, by doing so we do not 
lose generality, since we will be working in the interaction 
picture where only the frequencies Q n appear in relation 
to the state of the bath [Eq. ([12])]. Since the ratios of 
these frequencies and the temperature of the bath occur 
in the equations, only their values relative to the temper- 
ature are of interest. Therefore, henceforth we will omit 



the units of frequency and temperature and will treat 
these quantities as dimensionless. 

The interaction picture is defined as the transformation 
of any operator 

A i ^ A(t) = exp{iH a t)Aexp(-iH t), (4) 

where Hq = Hs + Hb ■ The interaction Hamiltonian Hi 
chosen here is invariant under this transformation since it 
commutes with Hq. [Note that in the next subsection, to 
simplify our calculations we redefine Hs and H'j (whence 
H'j becomes Hj), but this does not alter the present anal- 
ysis.] All the quantities discussed in the rest of this article 
are assumed to be in the interaction picture. 

The dynamics can be described using the superopera- 
tor notation for the Liouville operator 

Cp(t) = -i[Hi,p(t)], (5) 

where p(t) is the density matrix for the total system in 
the Hilbert space Hs ® Hb- The dynamics is governed 
by the von Neumann equation 

j t p(t) = aCp(t) (6) 

and the formal solution of this equation can be written 
as follows: 

p(t) = exp(a£t)p(0). (7) 

The state of the system is given by the reduced density 
operator 

Ps(t) = Tr s {p(i)}, (8) 

where Tr^ denotes a partial trace taken over the bath 
Hilbert space Hb- This can also be written in terms of 
the Bloch sphere vector 

fv x (t)\ 

v(t) = v v (t) = Tr{Z Ps (t)}, (9) 

where a = (a x , a v ,a z ) is the vector of Pauli matrices. In 
the basis of a z eigenstates this is equivalent to 

p s (t) = -(I + v-a) 

If l + v z (t) v x (t) - iv y (t)\ , , 
2 \v x (t)+iv y (t) l-v z (t) )■ [iyJ) 

We assume that the initial state is a product state, i.e., 

p(0) = ps(0)®p s , (11) 

and that the bath is initially in the Gibbs thermal state 
at a temperature T 

p B = exp(-H B /kT)/Tr[exp(-H B /kT)], (12) 

where k is the Boltzmann constant. Since pb commutes 
with the interaction Hamiltonian H], the bath state is 
stationary throughout the dynamics: = pb- Fi- 

nally, the bath spectral density function is defined as 
usual as 

J(O) =^| ffn | 2 5(f!-a0- (13) 
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B. Exact solution for the system-spin dynamics 

Wc first shift the system Hamiltonian in the following 
way: 



H s i * Hs + OI, 
9 = Tr{^9n<PB}. 



(14) 



As a consequence the interaction Hamiltonian is modified 
from Eq. {]]) to 



Hj Hi = ao z ® B, 



where 



D 



91 



B ■ 



(15) 



(16) 



This shift is performed because now Tib[Hi, p(0)] = 0, 
or equivalently 



Tr B {Bp B } = 0. 



(17) 



This property will simplify our calculations later when 
we consider approximation techniques in Sec. III. Now, 
we derive the exact solution for the reduced density op- 
erator ps corresponding to the system. We do this in two 
different ways. The Kraus operator sum representation 
is a standard description of the dynamics of a system 
initially decoupled from its environment and it will also 
be helpful in studying the coarse-graining approach to 
the quantum semigroup master equation. The second 
method is computationally more effective and is helpful 
in obtaining analytical expressions for N ^> 1. 



1. Exact Solution in the Kraus Representation 

In the Kraus representation the system state at any 
given time can be written as 



Ps 



(18) 



where (3 = 1/kT. Here 



N 1 

^ = Ea(-i)S 



(21) 



is the energy of each eigenstate |^), where I — Z x Z 2 ■ ■ ■ In is 
the binary expansion of the integer I, and the partition 
function is Z = J^i ex P( — f3Ei). Therefore, the Kraus 
operators become 



K, 



\i exp(—itaEi<7 z )5ij, 



(22) 



where 



JV 



E t = (i\B\i) = J29n(-iy n - Tr{J2 9n<Ps}, (23) 



and Xi = exp(—[3Ei)/Z. Substituting this expression for 
Kij into Eq. (fTS]) and writing the system state in the 
Bloch vector form given in Eq. fTQj) , we obtain 



where 



(24) 



(25) 



v x (t) = v x (0)C(t)-v y (0)S(t), 
v y {t) = v x (0)S(t) + v y (0)C(t), 
v z {t) - v z (0), 



C(t) = ^\ i cos2aE i t, 

i 

S(t) = ^2\ i sin2aE i t. 



The equations ([24]) are the exact solution to the system 
dynamics of the above spin bath model. We see that 
the evolution of the Bloch vector is a linear combination 
of rotations around the z axis. This evolution reflects 
the symmetry of the interaction Hamiltonian which is 
diagonal in the z basis. By inverting Eqs. (pH)) for ^(0) 
or v y (0), we see that the Kraus map is irreversible when 
C(t) 2 + S(t) 2 = 0. This will become important below, 
when we discuss the validity of the time-convolutionless 
approximation. 



where the Kraus operators satisfy J^.. K^Kij = Is (24J. 
These operators can be expressed easily in the eigenbasis 
of the initial state of the bath density operator as 



= \/Xi{j\exp(-iHit)\i}, 



(19) 



where the bath density operator at the initial time is 
Pb(0) — Y^i For the Gibbs thermal state chosen 

here, the eigenbasis is the iV-fold tensor product of the 
a z basis. In this basis 



Pb 



E 



exp(-/3£V 



(20) 



2. Alternative Exact Solution 

Another way to derive the exact solution which is com- 
putationally more useful is the following. Since all ct^ 
commute, the initial bath density matrix factors and can 
be written as 



JY 



PB = 



exp (-| 



kT u n) 



^Tr[exp(-^<)] 

N , N 



1 



(j +/?„<)= n 



(26) 
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where 



: i -Wf) 



(27) 



and — 1 < /3 n < 1. Using this, we obtain an expression 
for 6 defined in Eq. (H]) 



n— 1 m— 1 



= E.9nTr{i«+A-0} II + #»<)} 

n— 1 m^n 
N 

= Y.9nPn- (28) 



The evolution of the system density matrix in the inter- 
action picture is 



p s (t) = Tr B { e -^V(0)e^ 4 } 



(29) 



In terms of the system density matrix elements in the 
computational basis {|0), |1)} (which is an eigenbasis of 
<j z in Hi = aa z ® B) , we have 

U\Ps(t)\k) = (j\Tr B {e- iH 't 

N 

x Ps(0) (g) Pm e iHlt }\k) 

m— 1 

= TrB{e~' Q ^' crZ ' : ^ s * 

iV 

x 01ps(0)|fc)(g)p me +M<fc|CTZ|fc>B *}. 

m— 1 

Let us substitute (j|cr z |j) = (— l)- 7 and rewrite 

e -ia{j\<r"\j}Bt = e -ia(-l) 3 '(EiIiSl<-e-f> 



-i(-l) 3 a( S(C 7f-^/)t 



Thus the final expression for the system density matrix 
elements is 

(j\ps(t)\k) = (j\ P s(0)\k)e^ a9t 

N 

x IT i cof i(2e jk ag n t) - ij3 n sm(2e jk ag n t)] . 

Notice that eoo = en = 0, hence the diagonal matrix 
elements do not depend on time as before: 

(0| PS (t)|0) = (0|ps(0)|Q), 
(l\ Ps (t)\l) = (l|p s (0)|l). 

For the off-diagonal matrix elements eoi = 1, em = — 1, 
and the evolution is described by 

(o\ Ps (t)\i) = (oMo)|i>/(t), 

(l\ps(t)\0) = (l|p 5 (0)|0)rW, (30) 



where 



iV 



/(f) = e l2a8t Y[ [cos(2ag n t) - i{3 n sin(2a 5 „<)] . (31) 

n=l 

In terms of the Bloch vector components, this can be 
written in the form of Eq. (|2"4"]> . where 

c(t) = (/(*)+ /*(*))/2, 

S(t) = (f(t)-f*(t))/2i. (32) 



C. Limiting cases 

1. Short Times 

Consider the evolution for short times where at <C 1. 
Then 

JV 

I [ [cos(2ag„i) ± i(3 n sin(2ag n t)] 



- * 


JV / 

= nV 


1-(1- ffi)sm 2 (2ag n t) 


Since all the matrices are diagonal, they commute and 


n=l 
N 






we can collect the terms by qubits: 




- 2(1 ~f3l)(ag n t) 2 


] 


(3\ps(t)\k) = (j\ PS (0)\k) 


n=l 






N 


W 1-2 


n=l 


t 2 


m— 1 


« exp[- 


2(crf) 2 Q 2 ], 





(33) 



Let us denote (— l)- 7 — (— l) fe = 2ejA;. The trace can be 
easily computed to be 

Fill Tr{e- i2e ^(»»<- £ z )*i(J + /3„<)} 

II e i2e ^ a ^' [cos(2e jfe ag„f) - i/3„ sin(2e jk ag n t)] . 



where (see Appendix [A"| 



N 



Q 2 = Tr{B 2 PB } = Y J 9l^-Pl) 



2J(fi) 



-oo 1 -fcosh^) 



an. 



(34) 
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Note that for the above approximation to be valid, we 
need 2(at) 2 Q 2 < 1. The total phase of f{t) in Eq. ^3TJ) 
is 

N 

« 20crf + ^(-/3„2a 5 „i) 

n=l 

= 20ai-2a^ 5 „/?Jt = O, (35) 

where we have used Eq. (|28|) . Thus, the off-diagonal 
elements of the system density matrix become 

P f(t) « ^(0)e- 2 ^) 2 ^, 

p^) « ^°(0)e- 2 ^) 2 ^. (36) 

Finally, the dynamics of the Bloch vector components 
are: 

v x , y {t) « ^,,(0)e- 2 ( Qt ) 2< 3 2 , 

= «,(t). (37) 

This represents the well known behavior 25] of the evo- 
lution of an open quantum system in the Zeno regime. 
In this regime coherence does not decay exponentially 
but is initially flat, as is the case here due to the vanish- 
ing time derivative of p*g~(t) at t — 0. As we will see in 
Sec. Ill, the dynamics in the Born approximation (which 
is also the second order time-convolutionless approxima- 
tion) exactly matches the last result. 



while the temperature factors (3 n determine their modu- 
lation depths. If the codomain of spectral density is not 
continuous, i.e. it can be split into nonoverlapping inter- 
vals Gj, j — 1, J, then Eq. (f3"Tj) can be represented in 
the following form: 

f{t) = e aa0t P 1 {t)P 2 {t)...P J {t), (38) 

where 

p A t ) = II [cos(2ag n t)-i/3 n sm(2ag n t)]. (39) 

In this case, if Gj are separated by large enough gaps, 
the evolution rates of different Pj (t) can be significantly 
different. This is particularly noticeable if one Pj(t) un- 
dergoes partial recurrences while another Pj/ (t) slowly 
decays. 

For example, one can envision a situation with two 
intervals such that one term shows frequent partial re- 
currences that slowly decay with time, while the other 
term decays faster, but at times larger than the recur- 
rence time. The overall evolution then consists in a small 
number of fast partial recurrences. In an extreme case, 
when one g n is much larger then the others, this results in 
an infinite harmonic modulation of the decay with depth 
dependent on /?„, i.e., on temperature. 



4- Alternating signs 



2. Large N 

When N 1 and the values of g n are random, then 
the different terms in the product of Eq. (|3Tj) are smaller 
than 1 most of the time and have recurrences at differ- 
ent times. Therefore, we expect the function f(t) to be 
close to zero in magnitude for most of the time and full 
recurrences, if they exist, to be extremely rare. When g n 
are equal and so are fl n , then partial recurrences occur 
periodically, independently of N. Full recurrences occur 
with a period which grows at least as fast as N. This can 
be argued from Eq. (|24p by imposing the condition that 
the arguments of all the cosines and sines are simultane- 
ously equal to an integer multiple of 2ir. When J(f2) has 
a narrow high peak, e.g., one g n is much larger than the 
others, then the corresponding terms in the products in 
Eq. (j3"T]) oscillate faster than the rate at which the whole 
product decays. This is effectively a modulation of the 
decay. 



3. Discontinuous spectral density co-domain 

As can be seen from Eq. (f3"Tj) , the coupling constants 
g n determine the oscillation periods of the product terms, 



If the bath has the property that every bath qubit m 
has a pair —to with the same frequency fi_ m = f2 m , 
but opposite coupling constant g- m = —g m , the exact 
solution can be simplified. First, /3_ m = j3 m , and 9 = 0. 
Next, Eq. (|3"Tj) becomes 

N/2 

fit) = Yl [cos(2ag m t) - i/3 m sm(2ag m t)] x 

m— 1 

[cos(2ag- m t) - ifi- m sm(2ag- m t)] 

N/2 

= Y[ [cos 2 (2ag m t) + (3 2 m sin 2 (2ag m t)}. (40) 

m— 1 

This function is real, thus Eq. (|32p becomes C(t) = 
f(t),S(t) = 0, so that v x (t) = v x {0)f(t) and v y (t) = 
v y (0)f(t). The exact solution is then symmetric under 
the interchange v x <-> v y , a property shared by all the 
second order approximate solutions considered below, as 
well as the post-Markovian master equation. The limit- 
ing case Eq. (|3"3"j) remains unchanged, and since Q 2 de- 
pends on g 2 , but not g n , it and all second order ap- 
proximations also remain unchanged. In the special case 
1.9m I — 9i th e exact solution exhibits full recurrences with 
period T = tr/ag. 
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III. APPROXIMATION METHODS 

In this section we discuss the performance of differ- 
ent approximation methods developed in the open quan- 
tum systems literature 0, ■ The corresponding master 
equations for the system density matrix can be derived 
explicitly and since the model considered here is exactly 
solvable, we can compare the appoximations to the exact 
dynamics. We use the Bloch vector representation and 
since the z component has no dynamics, a fact which is 
reflected in all the master equations, we omit it from our 
comparisons. 



A. Born and Born-Markov approximations 

Both the Born and Born-Markov approximations are 
second order in the coupling strength a. 



1. Born approximation 

The Born approximation is equivalent to a truncation 
of the Nakajima-Zwanzig projection operator method at 
the second order, which is discussed in detail in Sec. Ill 
B. The Born approximation is given by the following 
integro-differential master equation: 

Ps(t) = - f Tr B {[H I (t),[H I {s),ps(s)®p B }}}ds. (41) 
Jo 

Since in our case the interaction Hamiltonian is time- 
independent, the integral becomes easy to solve. We ob- 
tain 

p s (t) = -2a 2 Q 2 f (p s (s) - a z p s (s)a z )ds, (42) 
Jo 

where Q 2 is the second order bath correlation function 
in Eq. (f34|) . Writing ps(t) in terms of Bloch vectors as 
(I + v ■ a)/2 [Eq. {ID}], we obtain the following integro- 
differential equations: 

v x , y (t) = -4a 2 Q 2 [ v x , y (s)ds. (43) 
Jo 

These equations can be solved by taking the Laplace 
transform of the variables. The equations become 



sV x ^ y (s) - v x<y (0) = ~4a 2 Q : 



(44) 



where V x ^ y (s) is the Laplace transform of v x , y (t). This 
gives 

v Xt y(0)s 



s 2 + 4Q 2 a 2 ' 



(45) 



which can be readily solved by taking the inverse Laplace 
transform. Doing so, we obtain the solution of the Born 
master equation for our model: 



Note that this solution is symmetric under the inter- 
change v x v y , but the exact dynamics in Eq. (|2~4"|) docs 
not have this symmetry. The exact dynamics respects the 



symmetry: v x 



and 



which is a symmetry 



of the Hamiltonian. This means that higher order cor- 
rections are required to break the symmetry v x <-> v y in 
order to approximate the exact solution more closely. 

One often makes the substitution v x ^ v (t) for v XiV (s) in 
Eq. {43} since the integro-differential equation obtained 
in other models may not be as easily solvable. This ap- 
proximation, which is valid for short times, yields 



v x , y (t) = -<±a 2 Q 2 tv x ^ y (t), 



which gives 



,(*) 



/ (0)exp(-2Q 2 a 2 t 2 



(47) 



(48) 



i.e., we recover Eq. ([37} . This is the same solution 
obtained in the second order approximation using the 
time-convolutionless (TCL) projection method discussed 
in Sec. Ill B. 



2. Born-Markov approximation 

In order to obtain the Born-Markov approximation, we 
use the following quantities Q[Ch.3]: 

R(lo) = ]T P Ei <t z Pe 2 , 

E2 — E\ —U) 

T{u) = a 2 / e^ s Q 2 ds, 
Jo 

H L = J2T(u)R{w?R{u), (49) 



where T(uo) — (T(uj) — T(uj)*)/2i, Ei is an eigenvalue 
of the system Hamiltonian ffg, and is the projector 
onto the eigenspace corresponding to this eigenvalue. In 
our case Hg is diagonal in the eigenbasis of <r z , and only 
uj = is relevant. This leads to R(0) = a z and T(0) — 
a 2 J °° Q 2 dt. Since T(0) is real, we have T(0) = 0. Hence 
the Lamb shift Hamiltonian Hr, = 0, and the Lindblad 
form of the Born-Markov approximation is 



Ps(t) = j{a z p S (T z - p s ), 



(50) 



v x , y (t) = v XtV (0) cos(2a\/Q 2 't). 



(46) 



where 7 = L(0) + L(0)* = 2a 2 Q 2 dt. But note that 
Q 2 = Tr b{B 2 pb} does not depend on time. This means 
that r and hence 7 are both infinite. Thus the Born- 
Markov approximation is not valid for this model and 
the main reason for this is the time independence of the 
bath correlation functions. The dynamics is inherently 
non-Mar kovian . 

A different approach to the derivation of a Marko- 
vian semigroup master equation was proposed in [23| . In 
this approach, a Lindblad equation is derived from the 
Kraus operator-sum representation by a coarse-graining 
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procedure denned in terms of a phenomenological coarse- 
graining time scale r. The general form of the equation 
is: 



dp{t) 

m 



-i[{Q)r,p{t)) 



1 M 

x (Xa, )r([A a ,p(t)Al} + [A aP (t),Al]), 



a, /3=1 



where the operators Ao — I and A a ,a — 1,...,M form 
an arbitrary fixed operator basis in which the Kraus op- 
erators (US]) can be expanded as 



M 



K L — ^ ' bi a A a . 



(51) 



a=0 



The quantities Xa,/3(t) and Q(t) are defined through 

x a A t ) = J2 bi <*® b U t )> 



M 



and 



(X)r = - fx 
T JO 



(s)ds. 



(52) 



(53) 



(54) 



For our problem we find 



dp{t) 
dt 



= -iu[az, p{t)] + j(azp(t)az - p(t)), (55) 



where 



and 



7 



2r 



(l-C(r)) 



(56) 



(57) 



with C(t) and S{t) defined in Eq. (QSJ). In order for 
this approximation to be justified, it is required that the 
coarse-graining time scale r be much larger than any 
characteristic time scale of the bath 23]. However, in 
our case the bath correlation time is infinite which, once 
again, shows the inapplicability of the Markovian approx- 
imation. This is further supported by the performance 
of the optimal solution that one can achieve by varying 
t, which is discussed in Sec. IV. There we numerically 
examine the average trace-distance between the solution 
to Eq. (|55|) and the exact solution as a function of r. The 
average is taken over a time T, which is greater than the 
decay time of the exact solution. We determine an opti- 
mal r for which the average trace distance is minimum 



and then determine the approximate solution. The solu- 
tion of Eq. ([55]) for a particular r in terms of the Bloch 
vector components is 

v x (t) = v x {0)C T (t)+v y (0)S T (t) 

v y (t) = v y (0)C T (t) - v x (0)S T (t), (58) 

where C T (t) = e"^ 7 ") 4 cos(<D(t)£) and S T (t) = 
e r ( T )* sin(a)(r)t). The average trace distance as a func- 
tion of t is given by, 

D (pcxact, PCG) = ^^IPexact - PCG I 
T 

= ^^)-c(t))H(sw-W 

t=Q 

x yjv x (0) 2 +v y (0) 2 , (59) 

where pco represents the coarse-grained solution and 
where \X\ = y/xTX. The results are presented in 
Sec. IV. Next we consider the Nakajima-Zwanzig (NZ) 
and the time-convolutionless (TCL) master equations for 
higher order approximations. 



B. NZ and TCL master equations 

Using projection operators one can obtain approximate 
non-Markovian master equations to higher orders in at. 
A projection is defined as follows, 



Vp = Tr B {p}<E>p B , 



(60) 



and serves to focus on the "relevant dynamics" (of the 
system) by removing the bath (a recent generalization is 
discussed in Ref . [26| ) . The choice of ps is somewhat ar- 
bitrary and can be taken to be ps(0) which significantly 
simplifies the calculations. Using the notation introduced 
in 03, define 



(S) = vsv 



(61) 



for any superoperator S. Thus (S n ) denote the mo- 
ments of the superoperator. Note that for the Liouvil- 
lian superoperator, (£) = by virtue of the fact that 
Ttb{Bpb{0)} = (see [!])• Since we assume that the 
initial state is a product state, both the NZ and TCL 
equations are homogeneous equations. The NZ master 
equation is an integro-differential equation with a mem- 
ory kernel jV(£, s) and is given by 



ps(t)®PB= Af{t,s)p s (s)® p B ds. (62) 
Jo 

The TCL master equation is a time-local equation given 
by 



Ps(t)®p B =JC{t) PS {t)® p B - 



(63) 
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When these equations are expanded in at and solved we 
obtain the higher order corrections. When the interaction 
Hamiltonian is time independent (as in our case), the 
above equations simplify to 



Af(t, s)ps(s) ® p B ds = 



OO 

£« 

n=l 



and 



(64) 



IC(t) 



— +n— i 

En 1 
(n-l)\ 

n=l y > 



(C ri 



(65) 



for the NZ and TCL equations, respectively, where the 
time-ordered integral operator I n {t,s) is defined as 



I n (t,s) 



dii / dt 2 - 

n Jo 



ds. 



(66) 



The definitions of the partial cumulants and the 

ordered cumulants (C) oc are given in Refs. Hl|- 
For our model we have 



<£} pc - (C) oc = 0, 



(67) 



and 



(O oc = (0-3</: 2 ) 2 



(68) 



Explicit expressions for these quantities are given in Ap- 
pendix[B] Substituting these into the NZ and TCL equa- 
tions (|6"4"|) and (f6"5)) . we obtain what we refer to below as 
the NZro and TCLn master equations, with n = 2,3,4. 
These approximate master equations are, respectively, 
second, third and fourth order in the coupling constant 
a, and they can be solved analytically. The second order 
solution of the NZ equation (NZ2) is exactly the Born ap- 
proximation and the solution is given in Eq. (|46[) . The 
third order NZ master equation is given by 

p s (t) = -2a 2 Q 2 l 2 (t, S )(p s (s)~a z ps(s)<T z ) 

+ i4a 3 Q 3 l 3 (t, S )(cr z ps(s)- p S {s)a z ), (69) 

and the fourth order is 

p s (t) = -2a 2 Q 2 l 2 (t, S )(p s (s) - a z p s (s)a z ) 
+ i4a 3 Q 3 l 3 (t, s)(a z p s (s) ~ ps(s)a z ) 
+ 8a 4 (Q 4 - QtW, s){p s (s) - a z p s (s)a z ). 

(70) 



These equations are equivalent to, respectively, 6th and 
8th order differential equations (with constant coeffi- 
cients) and are difficult to solve analytically. The results 
we present in the next section were therefore obtained 
numerically. 

The situation is simpler in the TCL approach. The 
second order TCL equation is given by 



Ps it) = -c?m B {[Hi,[Hi,p 8 (t)2> p B $)\\} 
= -2a 2 tQ 2 (p s (t)~a z ps(t)a z ), 



(71) 



whose solution is as given in Eq. (|4"5)) in terms of Bloch 
vector components. For TCL3 we find 

p s (t) = ~2a 2 tQ 2 (p s (t)~a z ps(t)a z ) 
t 2 

+ ^Q^ 3 -{a zPs {t)- p s (t)a z ), (72) 
and for TCL4 we find 



p s (t) = [-2a 2 tQ 2 + (8Q 4 -24Q^)a 4 -] 
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x (Ps(t) - cr z ps(t)a z ) 



+ 4iQ 3 a 3 -(<J z ps{t) - ps(t)cr z ). (73) 
These equation can be solved analytically, and the solu- 



(£ 2 ) 

\ 1 pc 


= (C 2 ) 


given by 


(c 2 ) 

\ / oc 


= (c 2 ) 


v x (t) = 


(C 3 ) 

\ 1 pc 


= <£»> 


Vy(t) = 




= <0 




<o PC = <o - 


-{C 2 f 


where g(t) 



f n (at) [v y {0) cos( ff (t)) - v x (0) sm(g(t))] 



(74) 



where g(t) = 4Q 3 a 3 t 3 /3, f 3 (at) = exp(-2Q 2 a 2 t 2 ) 
(TCL3) and f 4 {at) = exp(-2Q 2 a 2 t 2 + (2Q 4 - 
6Qj)aH 4 /3) (TCL4). It is interesting to note that the 
second order expansions of the TCL and NZ master equa- 
tions exhibit a v x <-> v y symmetry between the compo- 
nents of the Bloch vector, and only the third order cor- 
rection breaks this symmetry. Notice that the coefficient 
of a 3 does not vanish in this model unlike in the one con- 
sidered in because both (£ 3 ) pc ^ and (C 3 ) oc ^ 
and hence the third order (and other odd order) approx- 
imations exist. 



C. Post-Markovian (PM) master equation 

In this section we study the performance of the post- 
Markovian master equation recently proposed in [131 ]: 



dp® 
Bt 



= V 



dt'k(t') exp(Vt')p{t - t'). (75) 



This equation was constructed via an interpolation be- 
tween the exact dynamics and the dynamics in the 
Markovian limit. The operator T> is the dissipator in 
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the Lindblad equation (|50p . and k(t) is a phenomenolog- 
ical memory kernel which must be found by fitting to 
data or guessed on physical grounds. As was discussed 
earlier, the Markovian approximation fails for our model, 
nevertheless, one can use the form of the dissipator we 
obtained in Eq. ([50]) 



Vp = a z pa z 



(76) 



It is interesting to examine to what extent Eq. (|75[) can 
approximate the exact dynamics. As a measure of the 
performance of the post-Markovian equation, we will take 
the trace-distance between the exact solution p oxac t{t) 
and the solution to the post-Markovian equation pi(t). 
The general solution of Eq. (|75l) can be found by express- 
ing p(t) in the damping basis [29f and applying a Laplace 
transform fl3( ] . The solution is 

Pit) = E - E MLiP(t))Ri, (77) 



where 



Hi(t) = Lap 



s - X l k(s - \. t ) 



m(0) = Zi(t)in(p), 



(78) 

(Lap 1 is the inverse Laplace transform) with k being the 
Laplace transform of the kernel k, {Li} and {Ri} being 
the left and right eigenvectors of the superoperator T>, 
and Xi the corresponding eigenvalues. For our dissipator 
the damping basis is {LJ = {R { } = {^=, ^J=, 

and the eigenvalues are {0, —2, —2, 0}. Therefore, we can 
immediately write the formal solution in terms of the 
Bloch vector components: 



J x,y 



iV (t) = Lap 1 



1 



s + 2k(s + 2) 



w x ,„(0) = £(t)v XiV (0). 



(79) 



We see that v x (t) has no dependence on v y (0), and nei- 
ther does v v (t) on v x (0), in contrast to the exact solution. 
The difference comes from the fact that the dissipator D 
does not couple v x (t) and v v (t). This reveals an inher- 
ent limitation of the post-Markovian master equation: it 
inherits the symmetries of the Markovian dissipator T>, 
which may differ from those of the generator of the exact 
dynamics. In order to rigorously determine the optimal 
performance, we use the trace distance between the exact 
solution and a solution to the post-Markovian equation: 



D( Pcxact (t), Pl (t)) = -^( C {t)-at)? + s{ty 



x ^v x (o)z + v y (oy 



(80) 



Obviously this quantity reaches its minimum for = 
C(t),\/t independently of the initial conditions. The 



kernel for which the optimal performance of the post- 
Markovian master equation is achieved, can thus be for- 
mally expressed, using Eq. ([79|) . as: 



1 



'Lap 



1 



(81) 



Lap(C(t)) 

It should be noted that the condition for complete pos- 
itivity of the map generated by Eq. ([75|) . J2i£i(t)Lf <B> 
Ri > [H, amounts here to \£,(t)\ = \C(t)\ < 1, which 
holds for all t. Thus the minimum achievable trace- 
distance between the two solutions is given by 

(t), Pl (t)) = ±S(t)y/v x (0)*+v y (0)*. (82) 

The optimal fit is plotted in Sec. IV. 

Finding a simple analytical expression for the optimal 
kernel Eq. (|8ip seems difficult due to the complicated 
form of C(t). One way to approach this problem is to 
expand C(t) in powers of at and consider terms which 
give a valid approximation for small times at -C 1. For 
example, Eq. (|33[) yields the lowest non-trivial order as: 



C 2 (t) 



1 



2Q 2 a 2 t 2 



0(aH 



(83) 



Note that this solution violates the complete positivity 
condition for times larger than t — 1/ ' a\J2Q 2 . The cor- 
responding kernel is: 



k 2 (t) = 2a 2 Q 2 e 2t cos1i(2a 



(84) 



Alternatively we could try finding a kernel that matches 
some of the approximate solutions discussed so far. For 
example, it turns out that the kernel 



k-NZ 2 (t) 



2a 2 Q 2 e 2t 



(85) 



leads to an exact match of the NZ2 solution. Finding a 
kernel which gives a good description of the evolution of 
an open system is an important but in general, difficult 
question which remains open for further investigation. 
We note that this question was also taken up in the con- 
text of the PM in the recent study (30|, where the PM 
was applied to an exactly solvable model describing a 
qubit undergoing spontaneous emission and stimulated 
absorption. No attempt was made to optimize the mem- 
ory kernel and hence the agreement with the exact so- 
lution was not as impressive as might be possible with 
optimization. 



IV. COMPARISON OF THE ANALYTICAL 
SOLUTION AND THE DIFFERENT 
APPROXIMATION TECHNIQUES 

In the results shown below, all figures express the evo- 
lution in terms of the dimensionless parameter at (plot- 
ted on a logarithmic scale). We choose the initial condi- 
tion v x (0) = v y (0) — 1/V2 and plot only v x (t) since the 
structure of the equations for v x (t) and v y (t) is similar. 
In order to compare the different methods of approxima- 
tion, we consider various choices of parameter values in 
our model. 
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A. Exact Solution 

We first assume that the frequencies of the qubits in 
the bath are equal (fi„ = 1, Vn), and so are the coupling 
constants (g n = 1, Vn). In this regime, we consider large 
and small numbers of bath spins N = 100 and N = 4, 
and two different temperatures (3=1 and (3 = 10. Figs.[T] 
and [5] show the exact solution for N = 100 and N = 4 
spins, respectively, up to the second recurrence time. For 
each TV, we plot the exact solution for (3 = 1 and (3 = 10. 

We also consider the case where the frequencies fl n and 
the coupling constants g n can take different values. We 
generated uniformly distributed random values in the in- 
terval [—1,1] for both f2 n and g n . In Figs. and (0]) 
we plot the ensemble average of the solution over 50 ran- 
dom ensembles. The main difference from the solution 
with equal f2„ and g n is that the partial recurrences de- 
crease in size, especially as N increases. We attribute this 
damping partially to the fact that we look at the ensem- 
ble average, which amounts to averaging out the positive 
and negative oscillations that arise for different values of 
the parameters. The main reason, however, is that for a 
generic ensemble of random f2„ and g n the positive and 
negative oscillations in the sums (|25|) tend to average out. 
This is particularly true for large N, as reflected in Fig. [3] 
We looked at a few individual random cases for N = 100 
and recurrences were not present there. For N = 20 (not 
shown here), some small recurrences were still visible. 

We also looked at the case where one of the coupling 
constants, say gi, has a much larger magnitude than the 
other ones (which were made equal). The behavior was 
similar to that for a bath consisting of only a single spin. 
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FIG. 1: (Color online) Comparison of the exact solution at 
P = 1 and (3 = 10 for N = 100. 

In the following, we plot the solutions of different or- 
ders of the NZ, TCL and PM master equations and com- 
pare them for the same parameter values. 
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FIG. 2: (Color online) Comparison of the exact solution at 
(3 = 1 and (3 = 10 for N = 4. 
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FIG. 3: (Color online) Comparison of the exact solution at 
j3 = 1 and (3 = 10 for N = 100 for randomly generated g n 
and Q n . 



B. NZ 

In this subsection, we compare the solutions of differ- 
ent orders of the NZ master equation for fl n = g n = 1 . 
Fig. © shows the solutions to NZ2, NZ3, NZ4 and the 
exact solution for (3 = 1 and (3 = 10 up to the first re- 
currence time of the exact solution. For short times NZ4 
is the better approximation. It can be seen that while 
NZ2 and NZ3 are bounded, NZ4 leaves the Bloch sphere. 
But note that the approximations under which these so- 
lutions have been obtained are valid for at < 1. The 
NZ4 solution leaves the Bloch sphere in a regime where 
the approximation is not valid. For (3 = 10, NZ2 again 
has a periodic behavior (which is consistent with the so- 
lution), while the NZ3 and NZ4 solutions leave the Bloch 
sphere after small times. Fig. © shows the same graphs 
for N = 4. In this case both NZ3 and NZ4 leave the 
Bloch sphere for (3 = 1 and (3 = 10, while NZ2 has a 



11 



05- 
0.4- 



02- 
0.1 - 





10 10 



10" 10 10 10 



10 10 



FIG. 4: (Color online) Comparison of the exact solution at 
[3 — 1 and (3 — 10 for TV = 4 for randomly generated g n and 



periodic behavior. A clear conclusion from these plots 
is that the NZ approximation is truly a short-time one: 
it becomes completely unreliable for times longer than 
at < 1. 




FIG. 5: (Color online) Comparison of the exact solution, NZ2, 
NZ3 and NZ4 at f3 = 1 and /3 = 10 for TV = 100. The exact 
solution is the solid (blue) line, NZ2 is the dashed (green) 
line, NZ3 is the dot-dashed (red) line and NZ4 is the dotted 
(cyan) line. 



C. TCL 

Fig. plots the exact solution, TCL2, TCL3 and 
TCL4 at = 1 and j3 = 10 for TV = 100 spins and 
= 9n = 1. It can be seen that for (3 — 1, the TCL 
solution approximates the exact solution well even for 
long times. However, the TCL solution cannot repro- 
duce the recurrence behavior of the exact solution (also 
shown in the figure.) Fig. © shows the same graphs for 



FIG. 6: (Color online) Comparison of the exact solution, NZ2, 
NZ3 and NZ4 at = 1 and (3 = 10 for TV = 4. The exact 
solution is the solid (blue) line, NZ2 is the dashed (green) 
line, NZ3 is the dot-dashed (red) line and NZ4 is the dotted 
(cyan) line. 



TV = 4. In this case, while TCL2 and TCL3 decay, TCL4 
increases exponentially and leaves the Bloch sphere after 
a short time. This is because the exponent in the so- 
lution of TCL4 in Eq. (|74p is positive. Here again the 
approximations under which the solutions have been ob- 
tained are valid only for small time scales and the graphs 
demonstrate the complete breakdown of the perturbation 
expansion for large values of at. Moreover, the graphs re- 
veal the sensitivity of the approximation to temperature: 
the TCL fares much better at high temperatures. 

In order to determine the validity of the TCL approxi- 
mation, we look at the invertibility of the Kraus map de- 
rived in Eq. (JT5J) or equivalently Eq. (|2"5|) . As mentioned 
earlier, this map is non-invertible if C(t) 2 + S(t) 2 = for 
some t (or equivalently v x (t) — and v y {t) = 0). This 
will happen if and only if at least one of the (3 n is zero. 
This can occur when the bath density matrices of some 
of the bath spins are maximally mixed or in the limit of 
a very high bath temperature. Clearly, when the Kraus 
map is non-invertible, the TCL approach becomes invalid 
since it relies on the assumption that the information 
about the initial state is contained in the current state. 
This fact has also been observed for the spin-boson model 
with a damped Jaynes-Cummings Hamiltonian (2j. At 
the point where the Kraus map becomes non-invertible, 
the TCL solution deviates from the exact solution (see 
Fig. [9]). We verified that both v x and v y vanish at this 
point. 



D. NZ, TCL, and PM 

In this subsection, we compare the exact solution to 
TCL4, NZ4 and the solution of the optimal PM master 
equation. Fig. (fTU]) shows these solutions for TV = 100 
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FIG. 7: (Color online) Comparison of the exact solution, 
TCL2, TCL3 and TCL4 at (3 = 1 and (3 = 10 for N = 100. 
The exact solution is the solid (blue) line, TCL2 is the dashed 
(green) line, TCL3 is the dot-dashed (red) line and TCL4 is 
the dotted (cyan) line. Note that for = 1, the curves nearly 
coincide. 
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FIG. 8: (Color online) Comparison of the exact solution, 
TCL2, TCL3 and TCL4 at = 1 and = 10 for N = 4. 
The exact solution is the solid (blue) line, TCL2 is the dashed 
(green) line, TCL3 is the dot-dashed (red) line and TCL4 is 
the dotted (cyan) line. Note that for (3 = 1, TCL3, TCL4 and 
the exact solution nearly coincide. 



and (3 — 1 and (3 = 10 when £l n = g n = 1. Here we 
observe that while the short-time behavior of the exact 
solution is approximated well by all the approximations 
we consider, the long-time behavior is approximated well 
only by PM. 

For (3 = 1, NZ4 leaves the Bloch sphere after a short 
time while TCL4 decays with the exact solution. But 
as before, the TCL solution cannot reproduce the recur- 
rences seen in the exact solution. The optimal PM solu- 
tion, by contrast, is capable of reproducing both the de- 
cay and the recurrences. TCL4 and NZ4 leave the Bloch 
sphere after a short time for (3 = 10, while PM again 




FIG. 9: (Color online) Comparison of TCL2 and the exact 
solution to demonstrate the validity of the TCL approxima- 
tion for N — 4 and (3 = 1. The solid (blue) line denotes the 
exact solution and the dashed (green) line is TCL2. Note that 
the time axis here is on a linear scale. TCL2 breaks down at 
at ps 0.9, where it remains flat, while the exact solution has 
a recurrence. 



reproduces the recurrences in the exact solution. Fig. [TT] 
shows the corresponding graphs for N = 4 and it can be 
seen that again PM can outperform both TCL and NZ 
for long times. Figs. [T^] and Q]|] show the performance 
of TCL4, NZ4 and PM compared to the exact solution 
at a fixed time (for which the approximations are valid) 
for different temperatures ((3 € [0.01, 10]). It can be seen 
that both TCL4 and the optimal PM solution perform 
better than NZ4 at medium and high temperatures, with 
TCL4 outperforming PM at medium temperatures. The 
performance of NZ4 is enhanced at low temperatures, 
where it performs similarly to TCL4 (see also Figs. [llil 
and [TTj) . This can be understood from the short-time 
approximation to the exact solution given in Eq. (|37p . 
which up to the precision for which it was derived is also 
an approximation of NZ2 [Eq. (|4"rj|) ]. As discussed above, 
this approximation (which also coincides with TCL2) is 
valid when 2Q 2 {at) 2 < 1. As temperature decreases, 
so does the magnitude of Q2, which leads to a better ap- 
proximation at fixed at. Since NZ2 gives the lowest-order 
correction, this improvement is reflected in NZ4 as well. 

In Figs. [T4l and [l"5l we plot the averaged solutions over 
50 ensembles of random values for O n and g n in the inter- 
val [—1, 1]. We see that on average TCL4, NZ4 and the 
optimal PM solution behave similarly to the case when 
Q n = g n = 1. Due to the damping of the recurrences, es- 
pecially when N = 100, the TCL4 and the PM solutions 
match the exact solution closely for much longer times 
than in the deterministic case. Again, the PM solution 
is capable of qualitatively matching the behavior of the 
exact solution at long times. 
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FIG. 10: (Color online) Comparison of the exact solution, 
NZ4, TCL4 and PM at (3 = 1 and f3 = 10 for N = 100. 
The exact solution is the solid (blue) line, PM is the dashed 
(green) line, NZ4 is the dot-dashed (red) line and TCL4 is the 
dotted (cyan) line. Note that for = 1, TCL4, PM and the 
exact solution nearly coincide for short and medium times. 
Only PM captures the recurrences of the exact solution at 
long times. 
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FIG. 11: (Color online) Comparison of the exact solution, 
NZ4, TCL4 and PM at /3 = 1 and f3 = 10 for N = 4. The 
exact solution is the solid (blue) line, PM is the dashed (green) 
line, NZ4 is the dot-dashed (red) line and TCL4 is the dotted 
(cyan) line. Note that for f3 = 1, TCL4 and the exact solution 
nearly coincide for short and medium times. 

E. Coarse-graining approximation 

Finally, we examine the coarse-graining approximation 
discussed in Sec. III. We choose the time over which 
the average trace distance is calculated to be the time 
where the exact solution dies down. In Fig. [TH] we plot 
the coarse-grained solution for the value of r for which 
the trace distance to the exact solution is minimum. As 
can be seen, the coarse-graining approximation does not 
help since the Markovian assumption is not valid for this 



s 




qL 1 , , 1 1 1 1 ' ' ' 

10 10 10 10 

t 



FIG. 12: (Color online) Comparison of the exact solution, 
NZ4, TCL4 and PM at at = 0.1 for N = 100 for different 
/3 £ [0.01, 10]. The exact solution is the solid (blue) line, PM 
is the dashed (green) line, NZ4 is the dot-dashed (red) line 
and TCL4 is the dotted (cyan) line. 
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FIG. 13: (Color online) Comparison of the exact solution, 
NZ4, TCL4 and PM at at = 0.5 for N = 4 for different 
13 G [0.01, 10]. The exact solution is the solid (blue) line, PM 
is the dashed (green) line, NZ4 is the dot-dashed (red) line 
and TCL4 is the dotted (cyan) line. 



model. In deriving the coarse-graining approximation 
[lH one makes the assumption that the coarse-graining 
time scale is greater than any characteristic bath time 
scale. But the characteristic time scale of the bath is 
infinite in this case. 



V. SUMMARY AND CONCLUSIONS 

We studied the performance of various methods for ap- 
proximating the evolution of an Ising model of an open 
quantum system for a qubit system coupled to a bath 
bath consisting of N qubits. The high symmetry of the 
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FIG. 14: (Color online) Comparison of the exact solution, 
NZ4, TCL4 and PM at (3 = 1 and (3 = 10 for N = 100 for 
random values of g„ and Q n . The exact solution is the solid 
(blue) line, PM is the dashed (green) line, NZ4 is the dot- 
dashed (red) line and TCL4 is the dotted (cyan) line. Note 
that for (3 = 1 and (3 = 10, TCL4, PM and the exact solution 
nearly coincide. 
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FIG. 15: (Color online) Comparison of the exact solution, 
NZ4, TCL4 and PM at (3 = 1 and (3 = 10 for N = 4 for 
random values of g n and f2 n . The exact solution is the solid 
(blue) line, PM is the dashed (green) line, NZ4 is the dot- 
dashed (red) line and TCL4 is the dotted (cyan) line. Note 
that for = 1, TCL4, PM and the exact solution nearly 
coincide for short and medium times. 



model allowed us to derive the exact dynamics of the sys- 
tem as well as find analytical solutions for the different 
master equations. We saw that the Markovian approxi- 
mation fails for this model due to the time independence 
of the bath correlation functions. This is also reflected 
in the fact that the coarse-graining method [23| does 
not approximate the exact solution well. We discussed 
the performance of these solutions for various parame- 
ter regimes. Unlike other spin bath models discussed 
in literature (e.g., [l6]), the odd-order bath correlation 



FIG. 16: (Color online) Comparison of the exact solution 
and the optimal coarse- graining approximation for N = 50 
and (3=1. The exact solution is the solid (blue) line and 
the coarse-graining approximation is the dashed (green) line. 
Note the linear scale time axis. 



functions do not vanish, leading to the existence of odd- 
order terms in the solution of TCL and NZ equations. 
These terms describe the rotation around the z axis of 
the Bloch sphere, a fact which is reflected in the exact 
solution. We showed that up to fourth order TCL per- 
forms better than NZ at medium and high temperatures. 
For low temperatures we demonstrated an enhancement 
in the performance of NZ and showed that NZ and TCL 
perform equally well. We showed that the TCL approach 
breaks down for certain parameter choices and related 
this to the non-invertibility of the Kraus map describing 
the system dynamics. We also studied the performance 
of the post-Markovian master equation obtained in (l3j 
with an optimal memory kernel. We discussed possible 
ways of approximating the optimal kernel for short times 
and derived the kernel which leads to an exact fit to the 
NZ2 solution. It turns out that PM master equation per- 
forms as well as the TCL2 for a large number of spins and 
outperforms all orders of NZ and TCL considered here 
at long times, as it captures the recurrences of the exact 
solution. 

Our study reveals the limitations of some of the best 
known master equations available in the literature, in 
the context of a spin bath. In general, perturbative ap- 
proaches such as low-order NZ and TCL do well at short 
times (on a time scale set by the system-bath coupling 
constant) and fare very poorly at long times. These ap- 
proximations are also very sensitive to temperature and 
do better in the high temperature limit. The PM does 
not do as well as TCL4 at short times but has the distinct 
advantage of retaining a qualitatively correct character 
for long times. This conclusion depends heavily on the 
proper choice of the memory kernel; indeed, when the 
memory kernel is not optimally chosen the PM can yield 
solutions which are not as satisfactory [3fj| . 
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APPENDIX A: BATH CORRELATION FUNCTIONS 

Here we show how to calculate the bath correlation functions used in our simulations. The fc th order bath correlation 
function is defined as 

Q k = Tr{B k p B }, 

where B and ps were given in Eqs. p^]) and (fT2"|) . respectively. This yields: 



n I 

= j2 e ^^(i\(j: 9 n<-ei B n) 



i 



E &M z El \ l\C£9n< - ei B )\W\C£,9"<' ~ ei B)\m"\ ■ ■ ■ \l'"W"\(E9n>»<'» OlB) 



1,1',. ..,1"' 



E 6XP( f El) (Lan{i\<\i') - e)8 ll ,{Y J 9n'{i'K l \i") - e)6 n (E^»<z'"K»'IO - W*"i 



I n n' n'" 



or 



Q^I^^^oxrt), (Al) 



where Z = J2i ex p(—(3Ei) and the expressions for Ei and Ei were given in Eqs. ([2T|) and (|23|) . respectively. 

The above formulas are useful when the energy levels E\ and E\ are highly degenerate, which is the case for example 
when g n = g and £!„ = Q for all n. For a general choice of these parameters, it is computationally more efficient to 
consider 9 in the form l|2"8|) and the initial bath density matrix in the form (|2"rj)) . For example, the second order bath 
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correlation function is 

N N 

Q 2 = Tr{(]T g m (T z m - 9I)Q2 9n< - 01) PB } 

rn—l n—1 

N N 

= Tr { 9n9mK^PB}-29Tr(^2g n a^p B }- 

n,m=l n—1 



= Tr{ 9n 9m a z n a z m (g) ~(J + /?„<)} - 2 

n,m— 1 n—1 
N N 

= E Tr {5m ? «+/3 m /)}Tr{g„-K +/?„/)} } ( Tr{-(J + fitf)} + Tr{]T 9 2 iP b} - 2 



2 y n ' ' l " J-J- L 2 

N N N 

^2 9mPm9nPn ~ ^ 9nPn + ^ 5 « ~ ^ 
n,m— 1 n—1 n—1 



N 



= E^n(l -ft)- (A2) 

71=1 

Using the identity 1 — tanh 2 (— x/2) = 2/(1 + cosha;), this correlation function can be expressed in terms of the bath 
spectral density function [Eq. (fT3|) ] as follows: 

N 

Q2 = £ 



n=l 

°° n 

S(tl - n n )\g n \ 2 (l - tanh 2 (-— ))<M 
oo 2/eJ 

00 2J(Sl)cin 

_oo 1 + cosh(^)' 

Higher order correlation functions are computed analogously. 

APPENDIX B: CUMULANTS FOR THE NZ AND TCL MASTER EQUATIONS 

We calculate the explicit expressions for the cumulants appearing in Eq. (f6"5)) . needed to find the NZ and TCL 
perturbation expansions up to fourth order. 
Second order: 

(c 2 ) P = - , &B{[£r J ,[fl/,/j]]}®pB 

= - r lY B {Hjp-2H I pH I +pHj}®p B 

= -2Q 2 (PS - CF z p S (Tz) ® PB 

= p', (Bl) 

(C 2 ) 2 p = VC 2 VVC 2 Vp 
= V£ 2 Vp' 

= -2Q 2 {p' s - cr z p' s a z ) ® p B , 

where p' s = Tr B p' = -2Q 2 (p s - cr z pscr z ). Therefore 

(C 2 ) 2 p = -2Q 2 {(-2Q 2 (p s - (J z p s a z )) - a z (-2Q 2 (p s - a z p s a z ))a z } ® p B 

= 8Ql(ps — cr z ps<J z ) <g> pb- (B2) 
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Third order: 



Fourth order: 



(£> = il± B {[H I ,[H I ,[H I ,p]}]}®p B 

= iTr B {Hfp - ZHjpH! + 3H IP H] - pHj} ® p B 

= ±iQz{cFzPs- PsPz)® Pb- (B3) 



(£> = Tr s {[^ / ,[^ / ,[ J ff J ,[^ / ,p]]]]}®p B 

= Tr B {Hfp - 4H]pH! + QH 2 lP Hj - AHjpH* + pHJ} ® p B 

= 8Q 4 (ps -PzPsVz)® Pb- (B4) 
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